# load library
suppressPackageStartupMessages({
library(leaflet)
library(sf)
library(dplyr)
library(lwgeom)
library(purrr)
library(DT)
library(ggplot2)
library(mapview)
library(tmap)
library(ggspatial)
})
## Warning in fun(libname, pkgname): GEOS versions differ: lwgeom has 3.12.1 sf
## has 3.12.2
## Warning in fun(libname, pkgname): PROJ versions differ: lwgeom has 9.3.1 sf has
## 9.4.1
Data Source:
- NYC
Borough Boundaries
- NYC
Park Properties
- NYC
Street Centerline
target_crs <- 2263
# load data
borough_boundaries <- st_read("D:/GIS/R_compare/Borough Boundaries/geo_export_43cec398-2125-4dd7-bd74-d3b7e1dfe9f5.shp")
## Reading layer `geo_export_43cec398-2125-4dd7-bd74-d3b7e1dfe9f5' from data source `D:\GIS\R_compare\Borough Boundaries\geo_export_43cec398-2125-4dd7-bd74-d3b7e1dfe9f5.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 5 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -74.25559 ymin: 40.49613 xmax: -73.70001 ymax: 40.91553
## Geodetic CRS: WGS 84
parks <- st_read("D:/GIS/R_compare/Parks Properties_20241107/geo_export_d6f9bffe-7dce-4a9f-aae6-4ee8c4d82238.shp")
## Reading layer `geo_export_d6f9bffe-7dce-4a9f-aae6-4ee8c4d82238' from data source `D:\GIS\R_compare\Parks Properties_20241107\geo_export_d6f9bffe-7dce-4a9f-aae6-4ee8c4d82238.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 2051 features and 33 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -74.25609 ymin: 40.49449 xmax: -73.70905 ymax: 40.91133
## Geodetic CRS: WGS 84
streets <- st_read("D:/GIS/R_compare/NYC Street Centerline (CSCL)/geo_export_03186438-85ec-4a9d-9f79-f75e3ba6fa7e.shp")
## Reading layer `geo_export_03186438-85ec-4a9d-9f79-f75e3ba6fa7e' from data source `D:\GIS\R_compare\NYC Street Centerline (CSCL)\geo_export_03186438-85ec-4a9d-9f79-f75e3ba6fa7e.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 121949 features and 34 fields
## Geometry type: LINESTRING
## Dimension: XY
## Bounding box: xmin: -74.25496 ymin: 40.49788 xmax: -73.70002 ymax: 40.9151
## Geodetic CRS: WGS 84
streets <- streets[!st_is_empty(streets), ]
# Check and fix invalid geometries in parks
if (!all(st_is_valid(parks))) {
parks <- st_make_valid(parks) # Fix invalid geometries
}
Dynamic maps: Packages such as leaflet and mapview
in R support the creation of interactive maps that can be dynamically
displayed on web pages. In ArcGIS Pro, similar interactivity requires
the use of ArcGIS Online.
Flexible custom layers and layouts: ggplot2 and tmap
packages allow users to freely customize map layouts and layouts,
suitable for complex multi-layer drawing.
Leaflet is a powerful open-source JavaScript library used to create interactive maps. The R package ‘leaflet’ provides a convenient interface to integrate Leaflet’s capabilities within R, allowing users to build dynamic maps with ease. The package supports various map elements like tiles, markers, polygons, and popups, and is designed to be simple, efficient, and user-friendly.
# R: Dynamic Maps
# Create a Leaflet map with a basemap, legend, title, and scale
leaflet(parks) %>%
addTiles() %>% # Add online basemap
addPolygons(color = "green", weight = 1, fillOpacity = 0.5) %>%
addLegend(position = "bottomright", colors = "green", labels = "NYC parks", title = "Legend") %>%
addScaleBar(position = "bottomleft") %>%
addControl("<strong>NYC parks Map</strong>", position = "topright")
Advantages:
Interactive: The map generated by Leaflet can be panned and zoomed, supports click events and pop-up windows, and is very suitable for displaying interactive maps online.
Lightweight and fast: Leaflet is a lightweight library that loads
quickly and runs smoothly even on most browsers.
Open source and free: Compared with the commercial license fee of ArcGIS Pro, Leaflet is completely free and has rich community support and plug-in library.
Disadvantages:
Mapview is a fast visualization package in R, specifically for displaying geospatial data, based on Leaflet.
# Create an interactive map using mapview
map <- mapview(parks, legend = TRUE, layer.name = "NYC Parks", map.types = c("OpenStreetMap", "Esri.WorldImagery"))
# Show the map
map
Advantages:
Disadvantages:
ggplot2 is one of the most commonly used visualization libraries in R, suitable for making high-quality static maps.
# Create a static map using ggplot2
ggplot(borough_boundaries) +
geom_sf(aes(fill = boro_name), color = "black") + # 'boro_name' represents borough names
scale_fill_discrete(name = "Borough") + # Add legend for borough names
labs(
title = "NYC Borough Boundary",
caption = "Data Source: NYC OpenData"
) +
# Add scale bar and north arrow
annotation_scale(location = "bl", width_hint = 0.3) + # Add scale bar at the bottom left
annotation_north_arrow(location = "tl", which_north = "true", style = north_arrow_fancy_orienteering) + # Add north arrow at the top left
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
plot.caption = element_text(hjust = 0.5, size = 10),
legend.position = "bottom"
)
Advantages:
Disadvantages:
tmap is a package in R specifically for map making, suitable for static and interactive map production.
# Enable automatic fixing of invalid polygons (if needed)
tmap_options(check.and.fix = TRUE)
# Set tmap to static mode
tmap_mode("plot")
## tmap mode set to plotting
# Create a static map with a basemap
tm_shape(borough_boundaries) +
tm_polygons(title = "Borough Boundaries", col = "lightgrey", border.col = "black", legend.show = TRUE) +
tm_tiles("OpenStreetMap") + # Add a static basemap after all layers
tm_layout(
title = "NYC Map with Boroughs, Parks, and Streets",
title.size = 1.5,
title.position = c("center", "top"),
legend.position = c("left", "bottom"),
legend.bg.color = "white",
legend.frame = TRUE,
frame = FALSE,
bg.color = "lightblue",
inner.margins = c(0.05, 0.05, 0.05, 0.05)
) +
tm_scale_bar(position = c("left", "bottom")) +
tm_compass(type = "8star", position = c("right", "top")) +
tm_credits("Data Source: NYC Open Data", position = c("right", "bottom"))
# Set tmap to interactive mode
tmap_mode("view")
## tmap mode set to interactive viewing
# Create an interactive map
tm_shape(borough_boundaries) +
tm_polygons(title = "Borough Boundaries", col = "lightgrey", border.col = "black") +
tm_layout(
title = "NYC Map with Boroughs, Parks, and Streets",
title.size = 1.5,
title.position = c("center", "top"),
legend.position = c("left", "bottom"),
legend.bg.color = "white",
legend.frame = TRUE,
frame = FALSE,
bg.color = "lightblue",
inner.margins = c(0.05, 0.05, 0.05, 0.05)
) +
tm_basemap("OpenStreetMap") + # Add an online basemap
tm_scale_bar(position = c("left", "bottom"))
## legend.postion is used for plot mode. Use view.legend.position in tm_view to set the legend position in view mode.
Advantages:
Disadvantages:
Multi-layer clipping operations: R supports chaining operations, which allows you to clip multiple layers in sequence in a single code execution. In ArcGIS Pro, this requires multiple separate operations to complete.
Get the parks within Brooklyn.
# 1. Get the Brooklyn boundary
brooklyn_boundary <- borough_boundaries[borough_boundaries$boro_name == "Brooklyn", ]
brooklyn_boundary <- brooklyn_boundary[st_is_valid(brooklyn_boundary), ]
brooklyn_street<-streets[streets$borocode == 3, ]
brooklyn_parks<-parks[parks$borough == "B", ]
brooklyn_parks <- st_make_valid(brooklyn_parks)
brooklyn_street <- st_make_valid(brooklyn_street)
# 3. Perform multi-layer clipping to get parks near streets
parks_in_Brooklyn <- st_intersection(parks, brooklyn_boundary)
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
map1 <- mapview(brooklyn_parks, color = "blue", fill = "blue", alpha.regions = 1, lwd = 2, layer.name = "Brooklyn Parks by selection") +
mapview(parks_in_Brooklyn, color = "green", fill = "green", alpha.regions = 1, lwd = 2, layer.name = "Brooklyn Parks by clipping") +
mapview(brooklyn_boundary, color = "red", fill = NA, lwd = 2, layer.name = "Brooklyn Boundary")
map1
Flexible conditional expressions: R’s dplyr package can use more complex conditional expressions and functions to select data, supporting any custom conditions. In ArcGIS Pro, using the GUI to make such complex selections requires more steps and multiple conditional expressions.
Select areas in Brooklyn that have green space larger than a certain value.
brooklyn_boundary <- borough_boundaries %>% filter(boro_name == "Brooklyn")
# 2. parks in Brooklyn
suppressWarnings({parks_in_brooklyn <- st_intersection(parks, brooklyn_boundary)})
# 3. add area column
parks_in_brooklyn <- parks_in_brooklyn %>%
mutate(area_m2 = st_area(geometry))
# 4. select parks whose area > 1.5millions square meters
large_parks <- parks_in_brooklyn %>% filter(area_m2 > units::set_units(1500000, "m^2"))
# Plot the map and add a custom legend
ggplot() +
# Plot Brooklyn boundary (light grey)
geom_sf(data = brooklyn_boundary, aes(fill = "Brooklyn Boundary"), color = "lightgrey", size = 0.5, show.legend = "fill") +
# Plot large parks
geom_sf(data = large_parks, aes(fill = "large parks"), color = "green", size = 0.1, show.legend = "fill") +
# Add title and legend
labs(
title = "Parks in Brooklyn larger than 1.5 M (2024)",
caption = "Data Source: NYC Open Data",
fill = "Legend"
) +
# Set color mapping
scale_fill_manual(values = c(
"Brooklyn Boundary" = "lightgrey",
"large parks" = "green"
)) +
# Control border display in legend only
guides(fill = guide_legend(override.aes = list(color = "white", size = 0.5))) +
# Set theme style
theme_minimal(base_size = 14) +
theme(
plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
plot.subtitle = element_text(hjust = 0.5, size = 12),
plot.caption = element_text(hjust = 1, size = 10),
legend.position = "bottom"
)
Complex split operations: R supports splitting based on arbitrary geometry, and the st_split() function can split a layer based on custom shapes or other conditions. In ArcGIS Pro, splitting by arbitrary shapes is not always straightforward and requires complex layer editing. Compared with ArcGIS Pro, we can create a loop to split items in a list.
Use street data to segment Brooklyn park area.
large_parks <- large_parks[!st_is_empty(large_parks), ]
split_parks <- map(large_parks$geometry, ~ st_split(.x, streets))
split_parks_geom <- map(split_parks, ~ st_collection_extract(.x, "POLYGON"))
split_parks_combined <- do.call(c, split_parks_geom) %>%
st_sf()
map2 <- mapview(brooklyn_boundary,
col.regions = "lightgrey",
color = "darkgrey",
layer.name = "Brooklyn Boundary") +
mapview(split_parks_combined,
col.regions = "green",
color = "darkgreen",
layer.name = "Split Parks Combined")
## Warning in cbind(`Feature ID` = fid, mat): number of rows of result is not a
## multiple of vector length (arg 1)
map2
Join long streets’ attribute to the 2 largest park in Brooklyn
long_streets <- streets %>%
filter(st_length(geometry) > units::set_units(500, "m"))
joined_parks_streets <- st_join(large_parks, long_streets, join = st_intersects)
datatable(
joined_parks_streets,
options = list(pageLength = 5, autoWidth = TRUE, scrollX = TRUE),
caption = "Joined Parks and Streets Attributes Table"
)
Dynamic buffer: R supports the creation of dynamic buffers based on attributes, such as adjusting the buffer distance based on the attribute value of an object. In ArcGIS Pro, this attribute-based dynamic buffer needs to be implemented through a combination of attribute field settings and the buffer tool.
Find the buffer for the long streets based on their width.
buffered_polygon <- long_streets %>% st_buffer(dist = long_streets$st_width)
# Plot the map and add a custom legend
ggplot() +
# Plot Brooklyn boundary (light grey)
geom_sf(data = borough_boundaries, aes(fill = "NYC Boundary"), color = "darkgrey", size = 0.5, show.legend = "fill") +
# Plot buffer
geom_sf(data = buffered_polygon, aes(fill = "buffer zone"), color = "green", size = 0.1, show.legend = "fill") +
# Add title and legend
labs(
title = "Buffer for Long Street in NYC",
caption = "Data Source: NYC Open Data",
fill = "Legend"
) +
# Set color mapping for fill and border
scale_fill_manual(values = c(
"NYC Boundary" = "lightgrey",
"buffer zone" = "green"
)) +
guides(fill = guide_legend(override.aes = list(color = "darkgrey", size = 0.3))) +
# Set theme style
theme_minimal(base_size = 14) +
theme(
plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
plot.caption = element_text(hjust = 1, size = 10),
legend.position = "bottom"
)